% A Basic Mortensen-Pissarides Type of Search Model with SGU First-Order
% Local Approximation
% Written by: Orhan Torul
% Istanbul, Turkey, May 2017

clc;
clear all;
close all;

%%
%%
%Parameters

global beta delta xi kappa psi nu rhoz sigmaz sigma;
global zss;

beta=0.99;
delta=0.10; 
sigma=1;
xi=0.40;
kappa=0.50;
psi=0.25;
nu=0.40;
rhoz=0.95;
sigmaz=0.007;
%%
%%
%Steady-State
%c  : consumption
%u  : unemployment rate (n = 1-u)
%v  : vacancy postings
%q  : rate at which job vacations are filled
%p  : rate at which jobs are filled
%w  : wage
%R  : rate of return
%chi: unemployment benefit/home production determined as a ratio of w (0.4 or 0.95)
zss=0;

%y=[c; u ;v ;q ;p ;R ;w ;chi]
%y0=[0.255; 0.25; 0.5; 0.5; 0.5; 1.01; 0.25; 0.6];
y0=[0.87; 0.096 ; 0.277; 0.327 ; 0.945 ; 1.01 ; 0.9065 ; 0.3626];
%defining optimization options
options=optimset('Display','iter','MaxIter',10000, 'MaxFunEvals', 10000,'PlotFcns',@optimplotx,'tolfun',10e-25,'tolx',10e-25,'FunValCheck','on');
%calling for optimization
[out,fval,exitflag] = fsolve(@steady,y0,options); 

css=out(1);
uss=out(2);
vss=out(3);
thetass=out(3)/out(2);
qss=out(4);
pss=out(5);
Rss=out(6);
wss=out(7);
chiss=out(8);
%%
%Derivation of the Preliminaries of SGU Matrices

syms c u v p q chi w R z Xi;
syms c_1 u_1 v_1 q_1 p_1 z_1 w_1 Xi_1 R_1 chi_1; 

f1=-psi/(q*(1-delta))+Xi_1*((exp(z_1))-w_1+psi/q_1); 
%f1 is the job posting condition
f2=u*p-v*q;
%f2 is the equivalence of the match function definitions
f3=1-R*beta*((c_1/c)^(-sigma));
%f3 is the interest rate equation
f4=-w+nu*((exp(z))+psi*(v/u))+(1-nu)*chi;
%f4 is the wage as a result of wage bargaining
f5=-chi+0.40*w;
%f5 is the Shimer condition
f6=-(1-u_1)+(1-delta)*(1-u)+v*q;
%f6 is the evolution of the labor force 
f7=kappa*(u^xi)*(v^(1-xi))-v*q;
%f7 is the match function relation with q&v (and implicitly u&p)
f8=exp(z)*(1-u)+u*chi-c-psi*v;
%f8 is the aggregate resource constraint
f9=Xi_1-beta*((c_1/c)^(-sigma));
%f9 is the definition of discount rate Xi
f10=z_1-(1-rhoz)*zss-rhoz*z-randn(1,1)*sigmaz;
%f10 is the evolution of Solow residual/technology shock

F=[f1; f2; f3; f4; f5; f6; f7; f8; f9 ; f10];
%F is the vectorized form of all equilibrium conditions

y = [c; v; p; q; chi; w; R; Xi]; 
%y is the co-state vector as of time t
y_1= [c_1; v_1; p_1; q_1; chi_1; w_1; R_1; Xi_1];
%y_1 is the co-state vector as of time t+1
x = [u; z]; 
%x is the state vector as of time t
x_1= [u_1; z_1];
%x is the state vector as of time t+1


f_y=jacobian(F,y); 
f_y_1=jacobian(F,y_1);
f_x=jacobian(F,x);
f_x_1=jacobian(F,x_1);
%Above are the Jacobians to be used for SGU alghorithm.

%Setting the variables to their long-run (steady-state) values
c=css;
u=uss;
v=vss;
q=qss;
p=pss;
R=Rss;
w=wss;
chi=wss;
z=zss;

Xi=beta;
c_1=css; 
u_1=uss; 
v_1=vss;
q_1=qss;
p_1=pss;
z_1=zss;
w_1=wss;
chi_1=chiss;
R_1=Rss;
Xi_1=beta;

f_y=eval(f_y);
f_y_1=eval(f_y_1);
f_x=eval(f_x);
f_x_1=eval(f_x_1);
%Above are the steady-state values of the Jacobians

global a11g a12g a21g a22g a31g a32g a41g a42g a51g a52g a61g a62g a71g a72g a81g a82g;
global a11h a12h a21h a22h;
syms a11g a12g a21g a22g a31g a32g a41g a42g a51g a52g a61g a62g a71g a72g a81g a82g;
syms a11h a12h a21h a22h;

g=[a11g a12g; a21g a22g; a31g a32g; a41g a42g; a51g a52g; a61g a62g; a71g a72g; a81g a82g];
h=[a11h a12h; a21h a22h];


G =f_y_1*g*h + f_y*g + f_x_1*h + f_x;

Gvector=G(:);

Ginitial=[0.2;0.2;0.3;0.3;0.1;0.1;0.5;0.5;0.4;0.4;0.2;0.2;0.3;0.3;0.1;0.1;0.5;0.5;0.4;0.4;];
options=optimset('Display','iter','MaxIter',10000, 'MaxFunEvals', 10000,'PlotFcns',@optimplotx,'tolfun',10e-25,'tolx',10e-25,'FunValCheck','on');

[matrices,fval,exitflag] = fsolve(@(ppp) sgu(ppp,f_y_1,f_y,f_x_1,f_x),Ginitial,options); 

